#include <iostream>
#include <vector>
#include <fstream>
#include <cmath>
#include "Spline.hpp"
#include "Function.hpp"

const double PI = acos(-1.); // Define PI

class F1 : public Function {
public:
    double operator() (double x) const override {
        return 1 / (1 + 25 * x * x); // Example function
    }
};

class F2 : public Function {
public:
    double operator() (double x) const override {
        return std::sin(2 * PI * x); // Example function
    }
};

class F3 : public Function {
public:
    double operator() (double x) const override {
        return std::exp(-x * x); // Example function
    }
};

void testSpline(const Function& func, const std::string& funcName, std::vector<double> Knots, double step = 1e-2) {  
    // Generate linear piecewise polynomial objects using the function
    LinearPPForm Lin_PP(Knots, func);
    LinearBSpline Lin_B(Knots, func);
    // Output linear piecewise polynomial interpolation results to file and compare
    std::ofstream linearOutfile("output/T4_" + funcName + "_Linear.txt");
    std::ofstream EqOutfile("output/TestEquivalence.txt");
    for (double x = Knots.front(); x <= Knots.back(); x += step) {
        double y_PP = Lin_PP.evaluate(x);
        linearOutfile << x << " " << y_PP << std::endl;
    }
    linearOutfile << std::endl;
    for (double x = Knots.front(); x <= Knots.back(); x += step) {
        double y_B = Lin_B.evaluate(x);
        linearOutfile << x << " " << y_B << std::endl;
        EqOutfile << Lin_PP.evaluate(x) - y_B << std::endl;
    }
    linearOutfile.close();

    // Generate natural cubic spline objects using the function
    NaturalCubicPPForm Nat_PP(Knots, func);
    NaturalCubicBSpline Nat_B(Knots, func);
    // Output natural cubic spline interpolation results to file and compare
    std::ofstream naturalOutfile("output/T4_" + funcName + "_Natural.txt");
    for (double x = Knots.front(); x <= Knots.back(); x += step) {
        double y_PP = Nat_PP.evaluate(x);
        naturalOutfile << x << " " << y_PP << std::endl;
    }
    naturalOutfile << std::endl;
    for (double x = Knots.front(); x <= Knots.back(); x += step) {
        double y_B = Nat_B.evaluate(x);
        naturalOutfile << x << " " << y_B << std::endl;
        EqOutfile << Nat_PP.evaluate(x) - y_B << std::endl;
    }
    naturalOutfile.close();

    // Generate periodic cubic spline objects using the function
    PeriodicCubicPPForm Peri_PP(Knots, func);
    PeriodicCubicBSpline Peri_B(Knots, func);
    // Output periodic cubic spline interpolation results to file and compare
    std::ofstream periodicOutfile("output/T4_" + funcName + "_Periodic.txt");
    for (double x = Knots.front(); x <= Knots.back(); x += step) {
        double y_PP = Peri_PP.evaluate(x);
        periodicOutfile << x << " " << y_PP << std::endl;
    }
    periodicOutfile << std::endl;
    for (double x = Knots.front(); x <= Knots.back(); x += step) {
        double y_B = Peri_B.evaluate(x);
        periodicOutfile << x << " " << y_B << std::endl;
        EqOutfile << Peri_PP.evaluate(x) - y_B << std::endl;
    }
    periodicOutfile.close();

    // Generate complete cubic spline objects using the function
    double derivative_start = func.derivative(Knots.front());
    double derivative_end = func.derivative(Knots.back());
    CompleteCubicPPForm Comp_PP(Knots, func);
    CompleteCubicBSpline Comp_B(Knots, func);
    // Output complete cubic spline interpolation results to file and compare
    std::ofstream completeOutfile("output/T4_" + funcName + "_Complete.txt");
    for (double x = Knots.front(); x <= Knots.back(); x += step) {
        double y_PP = Comp_PP.evaluate(x);
        completeOutfile << x << " " << y_PP << std::endl;
    }
    completeOutfile << std::endl;
    for (double x = Knots.front(); x <= Knots.back(); x += step) {
        double y_B = Comp_B.evaluate(x);
        completeOutfile << x << " " << y_B << std::endl;
        EqOutfile << Comp_PP.evaluate(x) - y_B << std::endl;
    }
    completeOutfile.close();

    EqOutfile.close();
}

int main() {
    std::vector<double> evenKnots = {-1, -0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8, 1}; // Given uniform knots
    std::vector<double> UNevenKnots = {-1, -0.8, -0.7, -0.4, -0.2, 0, 0.3, 0.5, 0.6, 0.8, 1}; // Given non-uniform knots
    F1 func1;
    F2 func2;
    F3 func3;

    testSpline(func1, "F1", UNevenKnots);
    testSpline(func2, "F2", UNevenKnots);
    testSpline(func3, "F3", UNevenKnots);

    return 0;
}
